rmarkdown::render('main.Rmd', encoding = 'UTF-8', output_dir = "../docs")
rm(list = ls())
# Carregar os Pacotes
library(magrittr)
## Warning: package 'magrittr' was built under R version 3.4.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.4.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(glue)
## Warning: package 'glue' was built under R version 3.4.3
## 
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(lattice)
library(latticeExtra)
## Warning: package 'latticeExtra' was built under R version 3.4.1
## Loading required package: RColorBrewer
library(georob)
## Warning: package 'georob' was built under R version 3.4.2
## Loading required package: parallel
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.4.1
## 
## Attaching package: 'georob'
## The following object is masked from 'package:stats':
## 
##     step
library(sp)
library(mapview)
## Warning: package 'mapview' was built under R version 3.4.1
## Loading required package: leaflet
## Warning: package 'leaflet' was built under R version 3.4.1
library(raster)
## Warning: package 'raster' was built under R version 3.4.1
## 
## Attaching package: 'raster'
## The following object is masked from 'package:georob':
## 
##     cv
## The following object is masked from 'package:glue':
## 
##     trim
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:magrittr':
## 
##     extract
library(rmarkdown)
library(caret)
## Warning: package 'caret' was built under R version 3.4.1
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.4.1
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
## 
##     layer
# Sistemas de referência de coordenadas (Fonte: http://spatialreference.org/ref/epsg/)
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sirgas2000 <- sp::CRS('+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')

# Rampas de cores
col_soil_var <- topo.colors(100)

Caracterização dos dados

Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de \(Pinus taeda\) L. de 29 anos. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.

Na área de estudo foam identificados Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos, os quais representam a topossequência clássica da área de estudo (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).

Figura 1. Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).

Amostragem

Definimos uma malha amostral contendo 102 pontos para coleta de dados (Figura 1 b). Os pontos foram alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função \(clhs()\) implementado no pacote clhs. Foram consideradas como variáveis condicionantes da amostragem ELEV, VD, TWI, CNBL e DECLI as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).

Em cada ponto amostral foi medida a profundidade do solum (PF) e a altura das árvores (h). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras. Para representar h utilizamos o valor médio dos quatro indivíduos de \(Pinus taeda\) mais próximos de cada ponto amostral foram mensurados. Portanto, cada ponto de amostragem representa uma área (bloco) de aproximadamente 10 x 10 metros no campo.

pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)

sp :: coordinates(pontos) <- c('X' , 'Y')
sp :: proj4string(pontos) <- wgs84utm22s #coordenada referencia
pontos <- sp :: spTransform(pontos, wgs84utm22s) #transforma coordenada para wgs84

pontos@data
##     ID CLASSE AREIA_1 ARGILA_1    h   CAP  PF  PFd EA
## 1    1     RL     142      544 31.1 161.8  20  2.0 20
## 2    2     RL     159      432 32.8 151.3  10  1.0 20
## 3    3     RL      75      639 28.8 138.5  20  2.0 20
## 4    4     CX     121      628 29.5 162.3 100 10.0 30
## 5    5     RL     181      596 34.1 173.0  30  3.0 30
## 6    6     RL     121      647 32.9 146.8  30  3.0 30
## 7    7     RL      87      591 29.2 146.8  20  2.0 20
## 8    8     GX      93      564 35.2 159.3 100 10.0 30
## 9    9     RL     164      476 31.9 137.8  10  1.0 10
## 10  10     CX     135      558 33.1 139.3 100 10.0 50
## 11  12     RR     107      559 31.0 155.0  55  5.5 55
## 12  14     RL     178      539 29.6 136.8  10  1.0 10
## 13  15     RL     111      580 30.1 146.8  10  1.0 10
## 14  16     CX     100      578 34.5 148.3 100 10.0 40
## 15  17     RL     178      554 30.3 143.8  15  1.5 15
## 16  18     RL     141      573 24.0 132.8  10  1.0 10
## 17  19     RL     109      625 30.1 142.5  30  3.0 30
## 18  20     CX     137      375 34.0 155.0 100 10.0 35
## 19  21     CX     117      426 30.8 189.8 100 10.0 40
## 20  22     RL     128      470 31.7 151.8  20  2.0 20
## 21  23     CX     113      551 36.0 156.8 100 10.0 30
## 22  24     RL     134      534 31.5 140.8  15  1.5 15
## 23  25     CX     121      596 35.4 156.5 100 10.0 30
## 24  26     GX      99      423 32.9 166.0 100 10.0 40
## 25  27     RL     193      445 26.7 129.8  30  3.0 30
## 26  28     RR     137      550 34.7 147.3  50  5.0 50
## 27  29     RL     193      520 31.2 140.5  10  1.0 10
## 28  30     RR     156      551 31.3 139.0  40  4.0 40
## 29  31     CX      96      566 35.5 145.5 100 10.0 40
## 30  32     CX     151      534 34.0 146.3 100 10.0 30
## 31  33     CX     151      465 31.4 188.3  60  6.0 30
## 32  34     RL     152      504 25.5 148.5  10  1.0 10
## 33  35     CX     125      617 33.1 150.8 100 10.0 30
## 34  36     RR     119      627 28.5 142.8  30  3.0 30
## 35  37     CX     107      549 36.5 145.3 100 10.0 35
## 36  38     CX     157      555 34.2 139.3 100 10.0 30
## 37  39     CX     125      575 35.5 169.0  90  9.0 50
## 38  40     RL     147      498 29.0 149.0  15  1.5 15
## 39  41  RL/RR     107      580 32.6 149.6  30  3.0 30
## 40  42     CX     161      574 33.5 146.8 100 10.0 40
## 41  43     CX     127      614 34.1 142.8 100 10.0 30
## 42  44     RL     134      614 32.4 146.0  20  2.0 20
## 43  45     RR      73      330 29.1 127.8  40  4.0 40
## 44  46     CX      69      633 35.9 147.8 100 10.0 60
## 45  47     CX      81      664 33.7 136.5 100 10.0 40
## 46  48     RL     103      574 33.4 144.8  40  4.0 20
## 47  49     CX      90      643 32.1 146.0 100 10.0 40
## 48  51     RL     159      432 28.1 125.5  15  1.5 15
## 49  52     RL     181      596 30.5 123.8  35  3.5 20
## 50  53     RR     136      441 25.6 138.0  30  3.0 30
## 51  54     CX      93      641 30.0 148.5 100 10.0 30
## 52  56     RR      73      330 27.9 135.5  40  4.0 30
## 53  57     CX      40      614 36.2 171.5 100 10.0 40
## 54  58     RL     167      522 32.8 143.0  20  2.0 20
## 55  59     CX      59      432 28.8 137.8 100 10.0 40
## 56  60     RL     163      531 30.4 150.3  10  1.0 10
## 57  61     CX      51      453 32.5 176.5 100 10.0 50
## 58  62     RR      80      414 28.1 139.8  70  7.0 70
## 59  63     RL     202      512 31.5 147.0  30  3.0 20
## 60  64     CX     157      555 34.4 159.8 100 10.0 60
## 61  65     GX     127      630 30.6 152.3  70  7.0 30
## 62  66     GX      96      210 30.1 162.3  60  6.0 20
## 63  67     RR     100      471 33.1 154.0  50  5.0 50
## 64  68     RR     137      659 30.0 146.8  40  4.0 20
## 65  69     RL     229      495 24.4 140.0  10  1.0 10
## 66  70     RR     164      478 28.9 137.5  55  5.5 55
## 67  71     RR     120      633 30.0 170.3  60  6.0 20
## 68  72     CX      89      560 30.4 148.5  60  6.0 45
## 69  73     RL     152      563 27.2 160.0  10  1.0 10
## 70  74     RL     186      484 28.9 135.3  15  1.5 15
## 71  75     RL     149      549 32.2 130.0  30  3.0 30
## 72  76     RR      70      648 28.0 134.5  50  5.0 30
## 73  79     CX     163      482 30.8 151.3  60  6.0 20
## 74  80     CX     175      551 30.0 133.8  70  7.0 30
## 75  81     CX      58      499 28.9 175.0 100 10.0 40
## 76  82     CX      90      572 32.4 136.8 100 10.0 30
## 77  83     RL      75      639 28.0 130.0  20  2.0 20
## 78  84     RL     219      491 28.4 135.3  15  1.5 10
## 79  86     GX      62      425 28.8 154.5 100 10.0 40
## 80  87     CX      87      599 31.6 137.8  60  6.0 30
## 81  88     CX     117      533 31.6 132.0  70  7.0 30
## 82  89     CX      93      379 31.5 148.5  10  1.0 60
## 83  90     CX      68      600 31.5 137.3 100 10.0 40
## 84  91     RL      87      591 29.8 126.5  20  2.0 20
## 85  92     CX      92      576 28.5 160.0 100 10.0 40
## 86  93     CX      68      657 30.5 142.3  70  7.0 30
## 87  94     CX     109      611 29.8 139.5 100 10.0 40
## 88  95     CX      84      556 26.7 152.8  90  9.0 38
## 89  97     CX      84      556 31.0 140.3  80  8.0 40
## 90  99     OX      90      537 29.4 119.0 100 10.0 50
## 91 110     CX     100      558 30.8 153.3  60  6.0 20
## 92 114     RL     160      502 25.6 152.0  20  2.0 20
## 93 115     CX      80      549 29.1 134.8  90  9.0 60
## 94 122     CX      74      635 30.7 148.8  70  7.0 60
## 95 131     CX      30      714 29.2 146.8 100 10.0 40
## 96 132     CX     154      576 33.5 150.5 100 10.0 30

Verificar a normalidade dos dados

stats::shapiro.test(pontos$PFd)
## 
##  Shapiro-Wilk normality test
## 
## data:  pontos$PFd
## W = 0.83462, p-value = 5.627e-09
stats::shapiro.test(pontos$h)
## 
##  Shapiro-Wilk normality test
## 
## data:  pontos$h
## W = 0.98795, p-value = 0.5348

PF bimodal h normal

par(mfrow=c(1,2))

ex <- graphics::hist(pontos$PFd, breaks=10, xlab="", col="grey90", cex.axis=1,
          xlim=c(1,10),  ylim=c(0,37),ylab="",main=" ") 
    xfit<-base::seq(min(pontos$PFd),max(pontos$PFd), length=50) 
    yfit<-stats::dnorm(xfit,mean=mean(pontos$PFd),sd=sd(pontos$PFd))
   yfit<-  yfit*diff(ex$mids[1:2])*length(pontos$PFd) 
    lines(xfit, yfit, col="red", lwd=2) 
    mtext("Frequência",line=2.6, side=2, cex=1) 
    mtext("PF (dm)",line=2.6, side=1, cex=1) 
    graphics::rug((pontos$PFd), col="red")
     box()


ex2 <- graphics::hist(pontos$h, breaks=10, xlab="", col="grey90", cex.axis=1,
          xlim=c(23,38),  ylim=c(0,16),ylab="",main="") 
    xfit<-base::seq(min(pontos$h),max(pontos$h), length=50) 
    yfit<-stats::dnorm(xfit,mean=mean(pontos$h),sd=sd(pontos$h))
   yfit<-  yfit*diff(ex2$mids[1:2])*length(pontos$h) 
    lines(xfit, yfit, col="red", lwd=2) 
    mtext("Frequência",line=2.6, side=2, cex=1) 
    mtext("h (m)",line=2.6, side=1, cex=1) 
    graphics::rug((pontos$h), col="red")
     box()

Utilizei o MDE disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta \(reamostragem\) no software SAGA GIS.

A partir do MDE foram derivadas 12 variáveis topográficas utilizando a ferramenta \(Tarrain Analyses\) do software SAGA GIS.

ASP <- raster::raster("../data/Covars/ASP.tif")
CNBL <- raster::raster("../data/Covars/CNBL.tif")
DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
LS <- raster::raster("../data/Covars/LS.tif")
PLC <- raster::raster("../data/Covars/PLC.tif")
RSP <- raster::raster("../data/Covars/RSP.tif")
TPI <- raster::raster("../data/Covars/TPI.tif")
TRI <- raster::raster("../data/Covars/TRI.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
VD <- raster::raster("../data/Covars/VD.tif")
VDCN <- raster::raster("../data/Covars/VDCN.tif")
RFPF <- raster::raster("../data/Covars/predictionPF.tif")

A partir da função extract implementada no pacote raster extraí os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial “pontos”.

pontos$ASP <- raster::extract(ASP,  pontos)
pontos$CNBL <- raster::extract(CNBL, pontos)
pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$LS <- raster::extract(LS, pontos)
pontos$PLC <- raster::extract(PLC, pontos)
pontos$RSP <- raster::extract(RSP, pontos)
pontos$TPI <- raster::extract(TPI, pontos)
pontos$TRI <- raster::extract(TRI, pontos)
pontos$TWI <- raster::extract(TWI, pontos)
pontos$VD <- raster::extract(VD, pontos)
pontos$VDCN <- raster::extract(VDCN, pontos)
pontos$RFPF <- raster::extract(RFPF, pontos)

Além dos atributos do terreno, foi utilizado o mapa predito via random forest (RF) utilizando o comando \(rf()\) implementado no pacote caret, com número padrão de preditores a serem selecionados em cada nó padrão (mtry = padrão) e 1000 árvores \((ntrees = 1000)\), utilizando como preditoras as variáveis topográficas oriundas do MDE

sp::spplot(RFPF, scales = list(draw = TRUE))

Para fins de ordenamento de produção, a área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local.

Cada PIC é classificada em função da média de altura das 100 árvores de maior perímetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrimentos anuais é estabelecido o índices de sítio (IS). Esse índice é então atribuido à todo o talhão.

As parcelas da área foram classificação em 4 níveis, em que 1 corresponde ao sítio com melhor produtividade (Figura 2).

Figura 2 - Área de estudo, com indicação da área produtiva e parcelas de inventário contínuo contento os índices de sítio e respectivos valores de altura dominante

Os níveis de produtividade foram atribuidos aos polígonos e, naquelas em que não existe parcela de PIC, o valor de sítio foi estimado com base na altura das árvores amostradas á campo e das PCIs vizinhas.

pistola <- 
 raster::shapefile('../data/contorno/limite_projeto.shp',
                   stringsAsFactors = FALSE, encoding = 'UTF-8') %>%
sp::spTransform(wgs84utm22s)

Utilizei a função \(shapefile\) implementada no pacote raster para carregar o polígono da área de estudo - armazenado no objeto pistola e o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM).

ProdutividadePistola <- 
 raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
 sp::spTransform(wgs84utm22s)

sp::spplot(ProdutividadePistola, main = "Mapa de índices de Sítio")

#default.stringsAsFactors()
#str(ProdutividadePistola)

Usei a função \(over\) implementada no pacote sp para identificar o nível de produtividade do polígono dentro da qual cada observação se encontra e O resultado foi armazenado em uma coluna definida como Sitio no objeto mesmo objeto espacial pontos.

pontos$Sitio <- sp::over(x = pontos, y = ProdutividadePistola, na.omit = TRUE) %>% unlist()
pontos$Sitio <- as.factor(pontos@data$Sitio)
col_sitio <- terrain.colors(nlevels(ProdutividadePistola$ProdutividadePistola))

sp::spplot(
  ProdutividadePistola, scales = list(draw = T),
  main = 'Localização dos pontos') +
  lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
                  pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()

Parte I - Profundidade do solo

n = 96 observações para calibração + predição + validação (108 ha) Validação: validação cruzada (leave-one-out)

Modelo linear misto de variação espacial (georob)

Para utilizar esse modelo, supus que os dados são uma realização de um campo aleatório com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente. Assim, considerei como efeitos fixos da variação da PFd a predição da PF via floresta aleatória, atributos de terreno oriundos do MDE e os índices de sítio.

Ajuste do variograma:

limites <- seq(0, 2000, length.out = 15)
vario<- georob::sample.variogram(PFd ~ RFPF + ELEV + DECLI + Sitio,
    data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = c( "matheron"),
    mean.angle = TRUE, annotate.npairs = T) %>%
  plot()

lags <- seq(0, 2000, length.out = 15)
georob::sample.variogram(
 PFd ~ RFPF + ELEV + DECLI, data = pontos, locations = ~ X + Y, lag.dist.def = lags,
  xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180)) %>% 
    plot(type = "b", ylab = 'Semivariância', xlab = 'Distância de separação (km)')

Avaliada a evolução da semivariância nas diferentes direções, há evidência da existência de estruturas de autocorrelação espacial dependentes da direção, então o processo espacial é anisotrópico. Porém, assumi a isotropia do processo espacial e o semivariograma amostral gerado foi independente da direção.

Ajustei ao variograma amostral um modelo exponencial do variograma usando o método dos quadrados mínimos não-lineares ponderados, com ponderação definida conforme o método de Cressie. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = “BFGS”).

produzidos pelo otimizador a cada 10 iterações:

vario_fit <- 
  georob::fit.variogram.model(
  vario, variogram.model = 'RMexp', param = c(variance = 13, nugget = 10, scale = 70), 
  weighting.method = "cressie", method = "BFGS")
## 
## 
##                             variance       snugget        nugget         scale
##   Variogram parameters  0.000000e+00  0.000000e+00  0.000000e+00  3.204568e+22
## 
## 
##                             variance       snugget        nugget         scale
##   Variogram parameters  2.731246e+28  0.000000e+00  2.687536e-41  5.213753e+02
summary(vario_fit)
## 
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp", 
##     param = c(variance = 13, nugget = 10, scale = 70), weighting.method = "cressie", 
##     method = "BFGS")
## 
## Convergence in 212 function and 80 Jacobian/gradient evaluations
## 
## Residual Sum of Squares: 41.57611 
## 
## Residuals (epsilon):
##      Min       1Q   Median       3Q      Max 
## -2.37531 -1.66544 -0.61869  0.08033  1.99610 
## 
## Variogram:  RMexp 
##                 Estimate     Lower     Upper
## variance       1.312e+01 1.283e+01 1.341e+01
## snugget(fixed) 0.000e+00        NA        NA
## nugget         5.055e-04 1.088e-46 2.349e+39
## scale          6.635e+01 5.482e+01 8.029e+01
plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância')
lines(vario_fit, col = "red", lty = 'dashed')

Na geoestatística clássica, a variância não explicada é expressa em um único parâmetro nugget. Porém, no modelo linear misto podemos separar o parâmetro nugget em dois componentes: a variãncia devida aos erros de medida, modelada pelo parâmetro nugget e variância devida à variação espacial em pequena escala, modelada pelo parâmetro snugget. Considerando que durante a obteção dos dados de campo as informações de PF os valores foram arredondadas pelas em decímetros, considerei que a variância do erro de medida é igual a 0.5, ou seja, 62% do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget.

Feitas essa suposição, ajustei novamente o modelo novamente, agora mantendo ambos nugget e snugget fixos.

vario_fit_error <- georob::georob(
   PFd ~ RFPF, pontos, locations = ~ X + Y, variogram.model = 'RMexp', 
  param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']], 
            nugget = vario_fit$variogram.object[[1]]$param[['nugget']] * 0.5,
            snugget = vario_fit$variogram.object[[1]]$param[['nugget']] * 0.5,
            scale = vario_fit$variogram.object[[1]]$param[['scale']]),
  fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE),
  tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)
## 
## Call:georob::georob(formula = PFd ~ RFPF, data = pontos, locations = ~X + 
##     Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]], 
##     nugget = vario_fit$variogram.object[[1]]$param[["nugget"]] * 
##         0.5, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] * 
##         0.5, scale = vario_fit$variogram.object[[1]]$param[["scale"]]), 
##     fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE), 
##     tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
## 
## Tuning constant:  1000 
## 
## Convergence in 7 function and 6 Jacobian/gradient evaluations
## 
## Estimating equations (gradient)
## 
##                                   xi         scale
##   Gradient           : -1.073960e-02  1.324515e-02
## 
## Maximized restricted log-likelihood: -205.2834 
## 
## Predicted latent variable (B):
##      Min       1Q   Median       3Q      Max 
## -6.86420 -0.67666  0.09004  1.06093  6.73776 
## 
## Residuals (epsilon):
##        Min         1Q     Median         3Q        Max 
## -4.251e-04 -5.509e-05  7.842e-06  6.590e-05  4.148e-04 
## 
## Standardized residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4088 -0.4274  0.0635  0.5154  3.3177 
## 
## 
## Gaussian REML estimates
## 
## Variogram:  RMexp 
##                 Estimate     Lower  Upper
## variance       4.172e+00 3.103e+00  5.608
## snugget(fixed) 2.527e-04        NA     NA
## nugget(fixed)  2.527e-04        NA     NA
## scale          3.390e+01 1.829e+01 62.863
## 
## 
## Fixed effects coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.64442    0.73305  -4.972 2.98e-06 ***
## RFPF         0.15278    0.01117  13.674  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error (sqrt(nugget)): 0.0159 
## 
## Robustness weights: 
##  46 weights are ~= 1. The remaining 50 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.9423  0.9889  0.9956  0.9902  0.9978  0.9990
plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância', lty = 'dashed')
lines(vario_fit, col = "navyblue")
lines(vario_fit_error, col = "orange")

CRIAR O GRID: Suporte de predição pontual

grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
colnames(grid@coords) <- colnames(pontos@coords)
grid$Sitio <- sp::over(grid, ProdutividadePistola)
grid$DECLI <- raster::extract(DECLI, grid)
grid$ELEV <- raster::extract(ELEV, grid)
grid$RFPF <- raster::extract(RFPF, grid)


grid
## class       : SpatialPointsDataFrame 
## features    : 10004 
## extent      : 516813.9, 518852.5, 6908631, 6909820  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0 
## variables   : 4
## names       : Sitio,                DECLI,             ELEV,             RFPF 
## min values  :     1, 0.000142422635690309, 873.203247070312, 16.9653333333333 
## max values  :     4,    0.459370046854019, 936.466369628906, 94.5196666666667
grid <- 
  sp::SpatialPointsDataFrame(
    coords = grid@coords, 
    data = data.frame(grid),
    proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)

Argumento type: fazer a predição do sinal “signal” z, “response” se for y ou efeitos fixos “trend” se quiser predizer erro de medida, y = z + e, então response - signal = erro se tendência - predito = componente aleatório

Usar control.predict.georob(extended.output = TRUE)

pred_ponto <- predict(
  vario_fit_error, newdata = grid, type= "response", signif = 0.95,
  control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE
## Warning in points2grid(points, tolerance, round): grid has empty column/
## rows in dimension 1
str(pred_ponto)
## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
##   ..@ data       :'data.frame':  10004 obs. of  8 variables:
##   .. ..$ pred           : num [1:10004] 7.87 8.25 3.92 4.24 8.3 ...
##   .. ..$ se             : num [1:10004] 2.06 2.06 2.03 2.02 2.06 ...
##   .. ..$ lower          : num [1:10004] 3.8387 4.2183 -0.0577 0.2717 4.2687 ...
##   .. ..$ upper          : num [1:10004] 11.9 12.3 7.9 8.2 12.3 ...
##   .. ..$ trend          : num [1:10004] 7.94 8.31 4.08 4.4 8.4 ...
##   .. ..$ var.pred       : num [1:10004] 0.101 0.113 0.171 0.186 0.122 ...
##   .. ..$ cov.pred.target: num [1:10004] 0.0231 0.0263 0.1114 0.1332 0.0347 ...
##   .. ..$ var.target     : num [1:10004] 4.17 4.17 4.17 4.17 4.17 ...
##   .. ..- attr(*, "variogram.object")=List of 1
##   .. .. ..$ :List of 9
##   .. .. .. ..$ variogram.model: chr "RMexp"
##   .. .. .. ..$ param          : Named num [1:4] 4.17 2.53e-04 2.53e-04 3.39e+01
##   .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
##   .. .. .. ..$ fit.param      : Named logi [1:4] TRUE FALSE FALSE TRUE
##   .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
##   .. .. .. ..$ isotropic      : logi TRUE
##   .. .. .. ..$ aniso          : Named num [1:5] 1 1 90 90 0
##   .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
##   .. .. .. ..$ fit.aniso      : Named logi [1:5] FALSE FALSE FALSE FALSE FALSE
##   .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
##   .. .. .. ..$ sincos         :List of 6
##   .. .. .. .. ..$ co: num 6.12e-17
##   .. .. .. .. ..$ so: num 1
##   .. .. .. .. ..$ cp: num 6.12e-17
##   .. .. .. .. ..$ sp: num 1
##   .. .. .. .. ..$ cz: num 1
##   .. .. .. .. ..$ sz: num 0
##   .. .. .. ..$ rotmat         : num [1:2, 1:2] 1.00 -6.12e-17 6.12e-17 1.00
##   .. .. .. ..$ sclmat         : Named num [1:2] 1 1
##   .. .. .. .. ..- attr(*, "names")= chr [1:2] "" "f1"
##   .. ..- attr(*, "psi.func")= chr "logistic"
##   .. ..- attr(*, "tuning.psi")= num 1000
##   .. ..- attr(*, "type")= chr "response"
##   ..@ coords.nrs : num(0) 
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 516814 6908631
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   .. .. ..@ cellsize         : Named num [1:2] 8.68 8.68
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   .. .. ..@ cells.dim        : Named int [1:2] 236 138
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   ..@ grid.index : int [1:10004] 32416 32417 32139 32140 32180 32181 31903 31904 31905 31943 ...
##   ..@ coords     : num [1:10004, 1:2] 517534 517543 517178 517187 517534 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:10004] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "X" "Y"
##   ..@ bbox       : num [1:2, 1:2] 516810 6908627 518857 6909824
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "X" "Y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
spplot(pred_ponto)

at <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range()
sp::spplot(pred_ponto, zcol = c("lower", "pred", "upper"), main = "prediction")

sp::spplot(pred_ponto, zcol = 'se', main = "SE")

validacao <- georob::cv(vario_fit_error, nset = 95)
summary(validacao)
## 
## Statistics of cross-validation prediction errors
##       me      mede      rmse      made       qne      msse    medsse  
## 0.001978  0.129553  2.008528  1.379093  1.583218  1.072934  0.229126  
##     crps  
## 1.083096